library(survival)
library(KMsurv)
library(ggplot2)
library(ggpubr)
library(survminer)
library(plotly)
library(muhaz)
library(ggthemes)
data(burn)
burn1 <- burn
burn1 <-
data.frame(burn1, Treatment = factor(burn1$Z1, labels = c("Routine", "Cleansing")))
burn1 <-
data.frame(burn1, Gender = factor(burn1$Z2, labels = c("Male", "Female")))
burn1 <-
data.frame(burn1, Race = factor(burn1$Z3, labels = c("Nonwhite", "White")))
burn1 <- data.frame(burn1, PercentBurned = burn1$Z4)
burn1 <-
data.frame(burn1, SiteHead = factor(burn1$Z5, labels = c("NotBurned", "Burned")))
burn1 <-
data.frame(burn1, SiteButtock = factor(burn1$Z6, labels = c("NotBurned", "Burned")))
burn1 <-
data.frame(burn1, SiteTrunk = factor(burn1$Z7, labels = c("NotBurned", "Burned")))
burn1 <-
data.frame(burn1, SiteUpperLeg = factor(burn1$Z8, labels = c("NotBurned", "Burned")))
burn1 <-
data.frame(burn1, SiteLowerLeg = factor(burn1$Z9, labels = c("NotBurned", "Burned")))
burn1 <-
data.frame(burn1, SiteRespTract = factor(burn1$Z10, labels = c("NotBurned", "Burned")))
burn1 <-
data.frame(burn1, BurnType = factor(burn1$Z11, labels = c("Chemical", "Scald", "Electric", "Flame")))
burn1.surv <- with(burn1, Surv(T3, D3))
# plot(survfit(burn1.surv ~ Treatment, data = burn1),
# col = 1:2,
# lwd = 2)
# title("Time to Infection for Routine Care and Total Body Cleansing")
# legend(
# "topright",
# c("Routine Care", "Total Body Cleansing"),
# col = 1:2,
# lwd = 2
# )
#Lets observe the data
burn1
print(survdiff(burn1.surv ~ Treatment, data = burn1))
## Call:
## survdiff(formula = burn1.surv ~ Treatment, data = burn1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Treatment=Routine 70 28 21.4 2.07 3.79
## Treatment=Cleansing 84 20 26.6 1.66 3.79
##
## Chisq= 3.8 on 1 degrees of freedom, p= 0.05
# Plot the Kaplan Meier Curves
KMcurves <- survfit(burn1.surv ~ Treatment, data = burn1)
KMplot <-
ggsurvplot(KMcurves, ggtheme = theme_fivethirtyeight()) + labs(title = "KM Curves for Time to Infection for Routine Care \n and Total Body Cleansing")
ggplotly(KMplot[[1]])
cumHazPlot <-
ggsurvplot(
KMcurves,
fun = "cumhaz",
conf.int = FALSE,
palette = c("#5d8aa8", "#e32636"),
ggtheme =theme_fivethirtyeight()
) + ggtitle("Cumulative Hazard for Treatment Type")
#ggplotly(cumHazPlot[[1]])
#cumHazPlot
ggplotly(cumHazPlot$plot) # + geom_abline())
## Warning: plotly.js does not (yet) support horizontal legend items
## You can track progress here:
## https://github.com/plotly/plotly.js/issues/53
#Cloglog plot
cLogLogPlot <-
ggsurvplot(
KMcurves,
fun = "cloglog",
palette = c("#5d8aa8", "#e32636"),
ggtheme = theme_fivethirtyeight()
) + ggtitle("Complimentary Log-Log for Treatment")
ggplotly(cLogLogPlot[[1]])
## Warning: plotly.js does not (yet) support horizontal legend items
## You can track progress here:
## https://github.com/plotly/plotly.js/issues/53
#Building with time independent objects
survTimeIndep <- Surv(burn1$T3, burn1$D3)
#Build a cox model with just treatment
burn.cox.indep.T <- coxph(survTimeIndep ~ Treatment, data = burn1)
summary(burn.cox.indep.T)
## Call:
## coxph(formula = survTimeIndep ~ Treatment, data = burn1)
##
## n= 154, number of events= 48
##
## coef exp(coef) se(coef) z Pr(>|z|)
## TreatmentCleansing -0.5614 0.5704 0.2934 -1.914 0.0557 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## TreatmentCleansing 0.5704 1.753 0.321 1.014
##
## Concordance= 0.566 (se = 0.039 )
## Likelihood ratio test= 3.73 on 1 df, p=0.05
## Wald test = 3.66 on 1 df, p=0.06
## Score (logrank) test = 3.76 on 1 df, p=0.05
#Building with treatement and burnpercentage
# burn.cox.indep.TPB <- coxph(survTimeIndep ~ Treatment + PercentBurned, data = burn1)
# summary(burn.cox.indep.TPB)
#Building with treatement, burnpercentage, and burn Type
# burn.cox.indep.TPBbt <- coxph(survTimeIndep ~ Treatment + PercentBurned + BurnType, data = burn1)
# summary(burn.cox.indep.TPBbt)
#Building with treatement, burnpercentage, and burn Type + Respartory Tract
# burn.cox.indep.TPBbtResp <- coxph(survTimeIndep ~ Treatment + PercentBurned + BurnType + SiteRespTract, data = burn1)
# summary(burn.cox.indep.TPBbtResp)
#Building with treatement, burnpercentage, and burn Type + resp Tract and gender
# burn.cox.indep.TPBbtRespG <- coxph(survTimeIndep ~ Treatment + PercentBurned + BurnType + SiteRespTract + Gender, data = burn1)
# summary(burn.cox.indep.TPBbtRespG)
# Build a model with all of the above + race
summary(coxph(survTimeIndep ~ Treatment + PercentBurned + BurnType + SiteRespTract + Gender + Race, data = burn1))
## Call:
## coxph(formula = survTimeIndep ~ Treatment + PercentBurned + BurnType +
## SiteRespTract + Gender + Race, data = burn1)
##
## n= 154, number of events= 48
##
## coef exp(coef) se(coef) z Pr(>|z|)
## TreatmentCleansing -0.622888 0.536393 0.303140 -2.055 0.0399 *
## PercentBurned 0.002904 1.002909 0.007712 0.377 0.7065
## BurnTypeScald 1.582955 4.869324 1.099898 1.439 0.1501
## BurnTypeElectric 2.090905 8.092234 1.089444 1.919 0.0550 .
## BurnTypeFlame 0.915903 2.499032 1.024396 0.894 0.3713
## SiteRespTractBurned 0.216496 1.241718 0.359962 0.601 0.5475
## GenderFemale -0.559528 0.571479 0.402071 -1.392 0.1640
## RaceWhite 2.250113 9.488803 1.031151 2.182 0.0291 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## TreatmentCleansing 0.5364 1.8643 0.2961 0.9717
## PercentBurned 1.0029 0.9971 0.9879 1.0182
## BurnTypeScald 4.8693 0.2054 0.5639 42.0440
## BurnTypeElectric 8.0922 0.1236 0.9566 68.4549
## BurnTypeFlame 2.4990 0.4002 0.3356 18.6097
## SiteRespTractBurned 1.2417 0.8053 0.6132 2.5143
## GenderFemale 0.5715 1.7498 0.2599 1.2567
## RaceWhite 9.4888 0.1054 1.2575 71.6026
##
## Concordance= 0.738 (se = 0.039 )
## Likelihood ratio test= 24.8 on 8 df, p=0.002
## Wald test = 19.9 on 8 df, p=0.01
## Score (logrank) test = 23.44 on 8 df, p=0.003
# We can see from p-values that Treatment, Burn Type, and Race are the only variables with significance, so we will build a cox model using these
burn.cox <- coxph(survTimeIndep ~ Treatment + BurnType + Race, data = burn1)
summary(burn.cox)
## Call:
## coxph(formula = survTimeIndep ~ Treatment + BurnType + Race,
## data = burn1)
##
## n= 154, number of events= 48
##
## coef exp(coef) se(coef) z Pr(>|z|)
## TreatmentCleansing -0.6009 0.5483 0.2979 -2.018 0.0436 *
## BurnTypeScald 1.5785 4.8477 1.0865 1.453 0.1463
## BurnTypeElectric 2.1719 8.7745 1.0856 2.001 0.0454 *
## BurnTypeFlame 1.0035 2.7278 1.0163 0.987 0.3234
## RaceWhite 2.2842 9.8182 1.0260 2.226 0.0260 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## TreatmentCleansing 0.5483 1.8238 0.3058 0.983
## BurnTypeScald 4.8477 0.2063 0.5764 40.775
## BurnTypeElectric 8.7745 0.1140 1.0450 73.676
## BurnTypeFlame 2.7278 0.3666 0.3722 19.992
## RaceWhite 9.8182 0.1019 1.3144 73.340
##
## Concordance= 0.705 (se = 0.037 )
## Likelihood ratio test= 21.85 on 5 df, p=6e-04
## Wald test = 17.4 on 5 df, p=0.004
## Score (logrank) test = 20.61 on 5 df, p=0.001
LRT and Wald Test both provide evidence that there is a significant difference with patients with differrent treatments, burn types, and race.
# We check Residuals
#Schoenfield
burn.zph <- cox.zph(burn.cox)
ggcoxzph(
burn.zph[1],
se = FALSE,
font.main = 12,
ggtheme = theme_fivethirtyeight(),
point.col = "#5d8aa8",
)
ggcoxzph(
burn.zph[2],
se = FALSE,
font.main = 12,
ggtheme = theme_fivethirtyeight(),
point.col = "#5d8aa8",
)
ggcoxzph(
burn.zph[3],
se = FALSE,
font.main = 12,
ggtheme = theme_fivethirtyeight(),
point.col = "#5d8aa8"
)
# Martingale Residuals
burn.mart<- residuals(burn.cox, type = "martingale")
burn.lowess <- as.data.frame(lowess(burn1$T3, burn.mart))
# Plot Martingale
ggplotly(
ggplot() + aes(burn1$T3, burn.mart) + geom_point(col = "#5d8aa8") + labs(x = "Time", y = "Martingale Residuals", title = "Martingale Residuals vs. Time till Infection with Straphylocous Aureaus") + geom_line(data = burn.lowess, aes(x = x, y = y), col = "#388E3C") + theme_fivethirtyeight()
)
# We can see that we need a transformation, so we use log time
burn.lowess.log <- as.data.frame(lowess(log(burn1$T3), burn.mart))
# Plot Martingale
ggplotly(
ggplot() + aes(log(burn1$T3), burn.mart) + geom_point(col = "#5d8aa8") + labs(x = "Time", y = "Martingale Residuals", title = "Martingale Residuals vs. Time till Infection with Straphylocous Aureaus") + geom_line(data = burn.lowess.log, aes(x = x, y = y), col = "#388E3C") + theme_fivethirtyeight()
)
# Let us make time categorical
burn1["CatTime"] <-
as.factor(cut(burn1$T3,
c(0, 49, 100),
labels = c("<50", "50+")))
#Buld a new cox Model with new time variable
burn.cox.catTime <- coxph(survTimeIndep ~ Treatment + BurnType + Race + CatTime, data = burn1)
## Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
## Loglik converged before variable 6 ; coefficient may be infinite.
summary(burn.cox.catTime)
## Call:
## coxph(formula = survTimeIndep ~ Treatment + BurnType + Race +
## CatTime, data = burn1)
##
## n= 154, number of events= 48
##
## coef exp(coef) se(coef) z Pr(>|z|)
## TreatmentCleansing -6.226e-01 5.366e-01 2.966e-01 -2.099 0.0358 *
## BurnTypeScald 1.918e+00 6.809e+00 1.087e+00 1.765 0.0775 .
## BurnTypeElectric 2.200e+00 9.023e+00 1.086e+00 2.025 0.0428 *
## BurnTypeFlame 1.099e+00 3.000e+00 1.017e+00 1.080 0.2800
## RaceWhite 2.539e+00 1.266e+01 1.031e+00 2.462 0.0138 *
## CatTime50+ -1.906e+01 5.271e-09 3.880e+03 -0.005 0.9961
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## TreatmentCleansing 5.366e-01 1.864e+00 0.3000 0.9596
## BurnTypeScald 6.809e+00 1.469e-01 0.8095 57.2769
## BurnTypeElectric 9.023e+00 1.108e-01 1.0734 75.8500
## BurnTypeFlame 3.000e+00 3.333e-01 0.4088 22.0154
## RaceWhite 1.266e+01 7.897e-02 1.6775 95.5835
## CatTime50+ 5.271e-09 1.897e+08 0.0000 Inf
##
## Concordance= 0.743 (se = 0.035 )
## Likelihood ratio test= 37.73 on 6 df, p=1e-06
## Wald test = 18.65 on 6 df, p=0.005
## Score (logrank) test = 28.58 on 6 df, p=7e-05
### DOES NOT CONVERGE
# Now we add time dependent covariates
burn.cox.TD <- coxph(survTimeIndep ~ Treatment + BurnType + Race + T1 + D1+T2+D2 , data = burn1)
summary(burn.cox.TD)
## Call:
## coxph(formula = survTimeIndep ~ Treatment + BurnType + Race +
## T1 + D1 + T2 + D2, data = burn1)
##
## n= 154, number of events= 48
##
## coef exp(coef) se(coef) z Pr(>|z|)
## TreatmentCleansing -0.351516 0.703621 0.314807 -1.117 0.2642
## BurnTypeScald 1.440387 4.222331 1.097976 1.312 0.1896
## BurnTypeElectric 1.783797 5.952413 1.111258 1.605 0.1084
## BurnTypeFlame 0.831245 2.296176 1.021213 0.814 0.4157
## RaceWhite 2.304509 10.019254 1.029397 2.239 0.0252 *
## T1 0.001372 1.001373 0.021945 0.063 0.9502
## D1 -0.351020 0.703970 0.437191 -0.803 0.4220
## T2 0.001054 1.001054 0.015185 0.069 0.9447
## D2 -0.828456 0.436723 0.457322 -1.812 0.0701 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## TreatmentCleansing 0.7036 1.42122 0.3796 1.304
## BurnTypeScald 4.2223 0.23684 0.4909 36.320
## BurnTypeElectric 5.9524 0.16800 0.6742 52.553
## BurnTypeFlame 2.2962 0.43551 0.3103 16.993
## RaceWhite 10.0193 0.09981 1.3323 75.346
## T1 1.0014 0.99863 0.9592 1.045
## D1 0.7040 1.42052 0.2988 1.658
## T2 1.0011 0.99895 0.9717 1.031
## D2 0.4367 2.28978 0.1782 1.070
##
## Concordance= 0.754 (se = 0.036 )
## Likelihood ratio test= 30.51 on 9 df, p=4e-04
## Wald test = 25.62 on 9 df, p=0.002
## Score (logrank) test = 29.51 on 9 df, p=5e-04
# We remove the ones that are not significant to the model
# Now we add time dependent covariates
burn.cox.TD <- coxph(survTimeIndep ~ Treatment + BurnType + Race + D2 , data = burn1)
summary(burn.cox.TD)
## Call:
## coxph(formula = survTimeIndep ~ Treatment + BurnType + Race +
## D2, data = burn1)
##
## n= 154, number of events= 48
##
## coef exp(coef) se(coef) z Pr(>|z|)
## TreatmentCleansing -0.4218 0.6559 0.3034 -1.390 0.16455
## BurnTypeScald 1.5112 4.5324 1.0884 1.388 0.16500
## BurnTypeElectric 2.0049 7.4252 1.0867 1.845 0.06505 .
## BurnTypeFlame 0.8808 2.4129 1.0173 0.866 0.38658
## RaceWhite 2.3147 10.1214 1.0282 2.251 0.02438 *
## D2 -0.9012 0.4061 0.3386 -2.662 0.00777 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## TreatmentCleansing 0.6559 1.5246 0.3619 1.1888
## BurnTypeScald 4.5324 0.2206 0.5368 38.2649
## BurnTypeElectric 7.4252 0.1347 0.8825 62.4771
## BurnTypeFlame 2.4129 0.4144 0.3285 17.7206
## RaceWhite 10.1214 0.0988 1.3490 75.9383
## D2 0.4061 2.4625 0.2091 0.7885
##
## Concordance= 0.751 (se = 0.036 )
## Likelihood ratio test= 29.57 on 6 df, p=5e-05
## Wald test = 24.37 on 6 df, p=4e-04
## Score (logrank) test = 28.06 on 6 df, p=9e-05
drop1(burn.cox.TD, test = "Chisq")
# We check Residuals
#Schoenfield
burn.zph.TD <- cox.zph(burn.cox.TD)
ggcoxzph(
burn.zph.TD[1],
se = FALSE,
font.main = 12,
ggtheme = theme_fivethirtyeight(),
point.col = "#5d8aa8",
)
ggcoxzph(
burn.zph.TD[2],
se = FALSE,
font.main = 12,
ggtheme = theme_fivethirtyeight(),
point.col = "#5d8aa8",
)
ggcoxzph(
burn.zph.TD[3],
se = FALSE,
font.main = 12,
ggtheme = theme_fivethirtyeight(),
point.col = "#5d8aa8"
)
ggcoxzph(
burn.zph.TD[4],
se = FALSE,
font.main = 12,
ggtheme = theme_fivethirtyeight(),
point.col = "#5d8aa8"
)